Statistics of the work done by splitting a one-dimensional condensate 
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Motivated by experiments on splitting one-dimensional condensates, we study the statistics of 
the work done by a quantum quench in a bosonic system. We discuss the general features of the 
probability distribution of the work and focus on its behaviour at the lowest energy threshold, 
which develops an edge singularity. A formal connection between this probability distribution and 
the critical Casimir effect in thin classical films shows that certain features of the edge singularity are 
universal as the post-quench gap tends to zero. Our results are quantitatively illustrated by an exact 
calculation for non-interacting bosonic systems. The effects of finite system size, dimensionality, and 
non-zero initial temperature are discussed in detail. 
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I. INTRODUCTION 

The study of non-equilibrium dynamics in thermally 
isolated quantum system is one of the most active ar- 
eas of research in many-body quantum physics, mainly 
driven by ground-breaking experiments with ultra-cold 
atoms . Early experiments demonstrated the occur- 
rence of collapse and revival phenomena in systems of 
strongly correlated bosons, suggesting a high degree of 
quantum coherence despite the presence of strong many- 
body interactions . On the other hand, the peculiar form 
of the interactions and in particular their integrability 
has been conjectured to be the reason of the observed 
lack of thermalization in the non-equilibrium dynamics 
of quasi one-dimensional bosonic gases'. In turn, the 
closeness to integrability may cause the phenomenon of 
pre-thermalization t,-!l recently studied experimentally by 
looking at the dynamics of the statistics of the interfer- 
ence contrast of split quasi one-dimensional condensates 
on an atom chip . Atom chips, more than standard opti- 
cal lattices, are particularly suited for determining proba- 
bility distributions since they allow the separate observa- 
tion of the dynamics of each individual split condensate. 

Generically, the non-equilibrium dynamics is expected 
to depend crucially on the specific features of the sys- 
tem considered as well as on the way it is taken out of 
equilibrium . Among the many types of protocols stud- 
ied, the simplest one is perhaps the so-called quantum 
quench 1 ", consisting in a rapid change in the parame- 
ters of the hamiltonian of the quantum system. Not 
only can this protocol be rather easily implemented in 
experiments, but its theoretical description also reveals 
deep and interesting connections with boundary statisti- 
cal field theory and confined classical systems 11-1 . The 
effects of non-equilibrium protocols on many-body sys- 
tems can be characterized in a variety of ways, e.g., by 
studying the time dependence of correlation functions of 
local operators 2 . However, from a fundamental point of 
view, a quantum quench is analogous to a thermody- 



namic transformation and as such it should be character- 
ized by three macroscopic variables: the work W done on 
the system 14-16 , the entropy S produced 2,1 ' , and the heat 
Q possibly exchanged 1 ^. Focusing on the work W done 
on the system upon performing the quench, a common 
feature of both classical and quantum non-equilibrium 
processes is that W fluctuates among different realiza- 
tions of the same protocol 19-21 . Accordingly, its descrip- 
tion requires the introduction of a probability distribu- 
tion P(W). Studying P(W) offers several advantages 
since this distribution turns out to be connected to other 
quantities that have been extensively studied theoreti- 
cally in a variety of different physical contexts, such as 
the critical Casimir effect 1 - and the Loschmidt echo 1 '. 
These connections highlight emergent universal proper- 
ties of systems close to criticality 12,15 even for more gen- 
eral protocols 22,2 . On the other hand, work is a funda- 
mental observable and should be experimentally accessi- 
ble either directly by spectroscopic methods ' ' or in- 
directly by integrating the information extracted through 
the cloud expansion data of cold atom experiments. 

Our aim here is to discuss the qualitative features of 
P(W) expected for global quantum quenches of extended 
systems on the basis of the connections between quench 
dynamics and boundary statistical mechanics 10,12-14 . As 
we show below, the function P(W) consists of a super- 
position of peaks which correspond to the excited states 
and which gradually overlap to form a gaussian distribu- 
tion as the system size increases. While the distribution 
P(W) for "large" values of the work W depends on the 
microscopic details of the system, a universal behaviour 
might emerge 1 J in the form of a threshold edge singular- 
ity when W is comparable with the energy of the lowest 
excitations. This edge singularity, which has an exponen- 
tially small spectral weight in the thermodynamic limit, 
is a relevant feature in systems of finite size such as those 
currently investigated in experiments. Our general con- 
siderations are exemplified by detailed calculations of the 
statistics of the work in non-interacting bosonic systems 
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in d spatial dimensions subject to an abrupt change of 
the energy of the lowest excitation - which we will re- 
fer to as "mass" - from m to to. This simple model 
is of experimental relevance in various situations, e.g., 
for split quasi one-dimensional condensates in which the 
degrees of freedom associated with the relative phase of 
the two split components behave at low energies as non- 
interacting bosons. We also investigate the effect of a 
non-zero initial temperature on the properties of the 
statistics of the work. 

The rest of the presentation is organized as follows: In 
section II we introduce general concepts and discuss a 
possible experimental realization of the model we study. 
In section III we present the calculation of the statistics 
of the work for a quantum quench and derive some of 
its quantitative features. In section IV we explore in 
detail the connections with the critical Casimir effect. 
In section V we discuss the case of a quench originating 
from a thermal initial state. We summarize and discuss 
our results in section VI. 



II. EXPERIMENTAL MOTIVATION, MODEL 
AND QUALITATIVE RESULTS 

An interesting experimental realization of the quantum 
quench protocol, suitable for studying statistical fluctu- 
ations and probability distributions, is provided by the 
splitting of one-dimensional (lei) quasi-condensates, an 
experiment that has been performed extensively in the 
context of condensate interferometry over the last few 
years 8 . In these experiments, a system of ultra cold 
atoms trapped in a cigar-shaped potential and forming a 
Id quasi-condensate is split abruptly into two Id quasi- 
condensates, parallel to each other, by rapidly raising 
a potential barrier between them. The relative phase 
fluctuations of the two split condensates are the rele- 
vant low-energy degrees of freedom and their dynamics is 
described, under the circumstances discussed below, by 
a quadratic field theory. Indeed, while a Id system of 
bosons with short-range interactions may be described 
by the Lieb-Liniger model , its low-energy properties 
are effectively captured by the hydrodynamic or Lut- 
tinger liquid approach. The coupling of two such Id 
quasi-condensates due to tunnelling across the poten- 
tial barrier which separates them can in turn be mod- 
elled by a Josephson-like coupling. As a consequence, 
the hamiltonian H describing their relative phase (j){x) = 
(f>i(x) — </>2 (aO turns out to be 2 ' 1 ' 2, that of the sine-Gordon 
model 

H=^Jdx[U 2 + (d^) 2 - 4A cos(/^)] , (1) 

where H(x) is the momentum conjugate to the bosonic 
field 4>(x). The first two terms correspond to the Lut- 
tinger model while the last one — 4A cos(/3(/>) represents 
the Josephson coupling. The parameter A measures the 
strength of the coupling between the condensates. Indi- 



cating by M the mass of the bosons, by g their interac- 
tion strength (related to the scattering length) and by p 
their number density in three dimensions, for small val- 
ues of the dimensionless parameter 7 = Mg/(h 2 p) < 1, 
A is proportional to the hopping t between the two con- 
densates and to the density: A ~ t p. The parame- 
ter (3 in Eq. (1) is equal to \/27r /K where K is the 
Luttinger liquid parameter which is related solely to 7 
and varies from +00 for 7^1 (weakly interacting 
bosons) to 1 for 7 — > 00 (hard-core bosons or Tonks- 
Girardeau limit) . More specifically K increases like 
K ~ n/y/j(l — A /7/(27r)) -1 ' 2 upon decreasing 7 towards 
zero. In the limit of small interactions, i.e., /3 <C 1 (or, 
equivalently, large K or small 7), the system can be de- 
scribed by the quadratic approximation 

H = i J Ax [n 2 + {d x (t)) 2 + m 2 <j) 2 } + const. (2) 

characterized by a "mass" to = 2 A/3 2 for the quantum 
field 4>. Therefore we expect that the low-energy statis- 
tics of the work due to a partial, or total split of the 
two condensates can be described via this non-interacting 
bosonic model, which we will mostly focus on below. In 
d spatial dimensions, the hamiltonian in Eq. (2) is diag- 
onalizablc in independent momentum modes 



H(m) = 




where [<j>k,^k'\ = (2ir) d iS k ^. By J k we denote 
J BZ d d k/(2ir) d , i.e., the integral over the first Brillouin 
zone in the case of a lattice model, which becomes the 
whole d-dimensional space in the continuum limit. In the 
following we actually consider only this limit, although 
our discussion can be readily extended to the case on the 
lattice. Henceforth we assume a relativistic dispersion re- 
lation ojk(m) — \Jk 2 + to 2 and focus on quenches of the 
mass from mo to to. As we discussed above, this sim- 
ple model captures the physics of a number of physical 
systems, ranging from ideal harmonic chains to the low 
energy properties of interacting fcrmions and bosons in 
one dimension . 

Before describing qualitatively the main features of the 
statistics of the work done by quenching the mass in the 
quadratic field theory of Eqs. (2) and (3), we briefly recall 
some basic facts about the statistics of the work 14 ' 15 ' 19 ' 20 . 
In the case of a quench occurring from an initial state at 
zero temperature, the probability distribution P(W) of 
the work W done on the system by the quench is given 
by 

P(W) = J2*(W- E(n) + E (0))\(n\* o ) | 2 (4) 

n>0 

where \n) are the eigenstates of the post-quench hamil- 
tonian H with energies E(n) while Eo(0) is the energy 
of the ground state l^o) = |0)o of the pre-quench hamil- 
tonian Hq with eigenstates |n)o and eigenvalues Eo(n). 
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In the specific case of the quench we are interested in 
Hq = H(mo) while H = H{m) with H given in Eq. (3). 
Introducing the generating (or characteristic) function 20 



G(t) 



-iWt\ 



+00 



dW e 



-iWt 



P{W), (5) 



of the probability distribution P(W) one finds 

G(t) = (^ \e lHot e" lHt \^o) = e' i£ot (* |e"' iHt |*o)- (6) 

Performing a Wick rotation to imaginary "time" t t— > 
—iR, the amplitude ( v I / o|e _ ' f/fl i v I'o) can be interpreted as 
the partition function of the classical system correspond- 
ing to H, in a (d+l)-dimcnsional film of thickness R with 
boundary states described by l^o) on both boundaries of 
the film . Partition functions in such geometries are 
well-known in the study of the so-called critical Casimir 
effect 2 *' 2 '' and finite-size scaling ! "~ ! -'. In particular, the 
corresponding free energy can be decomposed into three 
contributions according to their scaling upon increasing 
R: 



InG(R) = -L d [f b R + 2f s + f c (R)} 



(7) 



where L — > 00 is the linear size of the system in 
all d dimensions parallel to the film. Here fb = 
liniL-s-oo AE gs /L d is the bulk free energy density, where 
AE gs = E(0) — E o (0); f s is the surface free energy as- 
sociated with each of the two (identical) boundaries of 
the film and fc{R) is the remaining finite-size contribu- 
tion which represents an effective interaction between the 
two boundaries and decays to zero for large R. Upon ap- 
proaching a critical point of the classical (and therefore 
quantum) system, this finite-size contribution acquires a 
universal character encoded in a scaling function of the 
ratio R/£ where £ is the bulk correlation length related 
to the inverse gap of the quantum system. This univer- 
sal contribution is known as the critical Casimir effect, 
which turned out to influence significantly the behaviour 
of, inter alia, soft matter systems and colloids at the mi- 
crometer and sub-micrometer scale - ^'" J ' \ In principle, 
G(t) in Eq. (6) can be calculated by analytic continuation 
of Eq. (7) to real "times" R = it, 



lnG(t) = -L d [if b t + 2f s + f c (it)}, 



(8) 



though particular care might be required in some cases 
for quenches across a quantum phase transition ' 1 . Ac- 
cordingly, when the post-quench hamiltonian is close to 
a critical point, G(t — > +00) — and therefore P(W) - 
is expected to acquire universal features which can be in- 
ferred from the critical Casimir effect. The work proba- 
bility distribution is then obtained via the inverse Fourier 
transform 



P(W) 



df 
2tt 



,+iWt 



G(t) 



(9) 



In the case of a quantum quench originating from a 
thermal initial state at a temperature /3 _1 , P(W) is given 
by a generalization of Eq. (4) : 

P(W) = J2 d ( W - E ( n )+ E °(no))Pn \(n\n ) \ 2 , (10) 



where the occupation probabilities p n of the pre-quench 
energy eigenstates are p n = e"^ E "' n '/2 n e"^ ^ and 
|n)o are the eigenstates of H with energies E (n). The 
corresponding generating function is therefore given by 

-PH iH t 



G{t) = 



Tr{c 



-iHt 



} 



Tr{e-^o} 



(11) 



We will refer to this quench originating from a thermal 
initial state as thermal quantum quench, in contrast to 
the case with zero initial temperature, which we will refer 
to as ordinary quantum quench. 




of G(t). 



FIG. 1. Probability distribution P(W) of the work W done 
on a single harmonic oscillator by a quench of its characteristic 
frequency from Wo = 1 to uj = 0.1. Vertical lines indicate 5- 
functions, the effective strength of which is represented in the 
vertical axis. 

The general features of the probability distribution 
P{W) associated with an ordinary quantum quench can 
be deduced already from Eqs. (4) and (8). It is instructive 
to consider first the case of a single quantum harmonic 
oscillator after a quench of its frequency from loq to uj. 
The associated P(W) consists of a sequence of equidis- 
tant Dirac <5-peaks corresponding to all the allowed tran- 
sitions from the initial ground state to the eigenstates 
of the post-quench hamiltonian, as shown in Fig. 1 for 
a specific choice of the values of ujq and lo. According 
to Eq. (4), the effective strength of each peak is set by 
the transition probability between the states which are 
involved in the transitions. Because of the behaviour of 
the wavefunctions of the harmonic oscillator under spa- 
tial inversion, only transitions to even excited states are 
allowed from the ground state (which is invariant under 
inversion) and consequently the spacing between consec- 
utive peaks is given by 2uj. The (5-peak at the lowest fre- 
quency corresponds to the transition to the new ground 
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FIG. 2. Probability density P(W) of the extensive work W done by an ordinary quantum quench in a system of Id non- 
interacting bosons with m = 1, mo = 0.1 and increasing system size L = 50, 75, 100, and 125. Sharp edge-singularity peaks 
appear at even integer multiples of m, superimposed on a wider bell-shaped distribution the mean and width of which increase 
upon increasing L. In particular, A indicates a contribution oc S(W) associated with the transitions between the pre- and 
post-quench ground states, whereas B indicates the lowest edge singularity, the shape of which is investigated here. 



state, since this requires the smallest possible amount of 
work. 

For the system of non-interacting (massive) bosons in 
Eq. (3) the transitions to excited states lead to the cre- 
ation of pairs of particles which, by conservation of mo- 
mentum, must have opposite momenta k and —k with en- 
ergies w_fc(m) = w&(m) > to. Therefore each of the pre- 
vious (5-peaks (except for the lowest one that corresponds 
to the transition from the pre- to the post-quench ground 
state) is now replaced by a continuous distribution with 
a lower threshold corresponding to the minimum energy 
cost for the creation of n pairs of particles with mass to, 
i.e., 2nm. This is clearly visible in Fig. 2, where we report 
P(W) for a given quench of the mass to of the model (3) 
in one spatial dimension, for various value of the size 
L of the chain. Given that P(W) vanishes identically 
below the extensive threshold AE gs ~ fbL d , in Fig. 2 
and in what follows we refer W to this value, so that 
the transition between the ground states of the pre- and 
post-quench hamiltonian give rise to a <5-peak located at 
W — 0. In passing we observe that AE gs is nothing but 
the work W lev that one would do at zero temperature 
= on the system during a transformation which 
changes the parameter to from its pre- to its post-quench 
value and which is reversible in the sense of thermody- 
namics. In fact, this reversible work W Tev is equal to the 
difference AF = Fp(m) — Fp(mo), where Fp(m) is the 



free energy of the system in contact with a thermal bath 
at temperature /3 _1 and hamiltonian H(m), which ren- 
ders the ground-state energy for ft — > oo. Accordingly, 
the mean of the quantity W — AE gs — W ~ W ICV stud- 
ied below is nothing but the so-called irreversible work 
Wirr' ! 5 (see, e.g., Ref. 17 where W lvr and the associated ir- 
reversible entropy production AS'i rr is studied for a ther- 
mal quench of the Ising chain) . 

In an interacting bosonic system, the pattern described 
above may no longer consist of equidistant peaks and 
it may involve excitations that are not simply pairs of 
particles with opposite momenta (as these properties are 
specific to non-interacting systems) but some of the fea- 
tures of P(W) are actually more robust. In particular, 
the low-energy part of P(W) always consists of a <5-peak 
at W = 0, with an amplitude given by the square J- of 
the so-called fidelity |(0|\&o)| between the pre- and post- 
quench ground state which, from Eq. (8), is 

J-=|(0|* }| 2 =e- 2L< ^. (12) 

In addition, provided the system has a gap to 7^ 0, a 
lower threshold emerges at W — 2m which reflects the 
fact that the lowest allowed excitations consist of pairs of 
quasiparticles with opposite momenta, each of which has 
energy at least equal to the energy gap to. Just above this 
threshold, P(W) turns out to exhibit an edge singularity 
whose profile is determined by the way fc in Eq. (7) de- 
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cays for R — > +00. Depending on the space dimensional- 
ity d of the system, this edge singularity (already present 
in the non-interacting case) displays different qualitative 
features, as shown in Fig. 3 for an instance of quench 
of the model in Eq. (3) and d = 1, 2, and 3 (from top 
to bottom). Further below we rationalize these features 
in full generality, by exploiting the connection with the 
critical Casimir effect. 

The fine structure of overlapping peaks (described 
above for the non-interacting case) fades out, however, 
upon increasing the system size L and P(W) gradually 
develops a highly peaked gaussian distribution with an 
extensive mean value oc L d and a relative width cx L~ d / 2 
which vanishes as L —¥ +00. As we show rigorously fur- 
ther below, this is a direct consequence of the fact that 
lnG(i) is extensive, i.e., proportional to L d , as can be 
seen from Eq. (8). The emergence of a gaussian distribu- 
tion as L increases is shown in Fig. 4 where, for the same 
system and quench as in Fig. 2, we report the distribu- 
tion P(w) of the intensive work w = W/L d for increasing 
values of L. While typical fluctuations around the mean 
(and most probable) value of w are indeed characterized 
by a gaussian distribution as L —¥ 00, large fluctuations 
display non-trivial and universal features 1 . 

For a thermal quantum quench, thermal excitations in 
all energy levels pre-exist in the initial state and they 
do not form only pairs of opposite momenta as in the 
case of the ordinary quench discussed above. Thus addi- 
tional peaks appear in P(W) corresponding to transitions 
from all the energy levels of the pre-quench hamiltonian 
to all those of the post-quench one, with an amplitude 
which decays exponentially upon increasing the corre- 
sponding energy difference. These peaks of thermal ori- 
gin are present also below the lower threshold — which 
is a characteristic feature of the quench originating from 
a zero-temperature initial state (i.e., the ground state, 
or more generally, a generic eigenstate of the pre-quench 
hamiltonian) — and they obscure the zero-temperature 
structure of P(W). In addition, for a thermal quench, 
the relation of P(W) with the statistical mechanics of a 
classical system in a film encoded in Eq. (8) is no longer 
valid. 



III. ORDINARY QUANTUM QUENCH 

In this section we show explicitly how the qualitative 
features mentioned above emerge in the case of the non- 
interacting bosonic model of Eq. (3) and we discuss in 
detail the scaling properties of P(W) as a function of 
the system size L. The calculation of G(t) for quenches 
of the mass in Eq. (3) is clearly simplified by the non- 
interacting nature of the problem. The initial state l^o) 
can be decomposed in the basis of eigenstates of the post- 
quench hamiltonian H (generated by the action of the 
post-quench creation operators a\ on the vacuum |0) of 
H) as a set of "Cooper pair" excitations; more specifically 
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FIG. 3. Lowest edge singularity of P(W) for a quench in 
the non-interacting bosonic system in Eq. (3) with m = 1 
and mo = 0.8 and, from top to bottom, spatial dimensional- 
ity d = 1, 2 and 3. The solid (blue) lines are the numerical 
evaluation of Eq. (9) with Eq. (16), while the dashed (red) 
lines correspond to the analytic expressions in, c.f., Eqs. (33) 
and (34). Here W is measured in units of m and, for con- 
venience, different values of L have been used in each plot. 
Note that, depending on the dimensionality d, P(W) diverges 
(d < 2), tends to a finite value (d — 2) or vanishes (d > 2) 
upon approaching the threshold W = 2m. The agreement 
between the numerical and the analytical curves is excellent 
for W < Am, where the analytic expression is indeed expected 
to reproduce the numerical data, as discussed further below 
in section IV. The noise in the blue curves, which becomes 
increasingly relevant as the dimensionality increases, is an ar- 
tifact of the numerical integration. 
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(see appendix A): 

l^o) = 11(1 - Ag) 1 / 4 exp 



*k f 



4)|0), (13) 



where 



A, 



wot + l^fe 



(14) 



Here and in what follows, u ak = ^(toq) and u> k = u> k (m) 
indicate the dispersion relations of the pre- and post- 
quench hamiltonian, respectively. This very special type 
of state is a so-called squeezed vacuum state and consists 
solely of pairs of particles of opposite momenta. As dis- 
cussed in appendix A, it emerges because the pre- and 
post-quench creation/annihilation operators are related 
by a linear Bogoliubov transformation. These transfor- 
mations and the resulting squeezed states characterize 
a large number of quantum quenches in non-interacting 
models, both bosonic and fermionic, while in interact- 
ing systems (that cannot be equivalcntly described by 
non-interacting ones) they are not of such generic util- 
ity but they emerge, instead, only in special cases of 
quenches ' . 

Due to the non-interacting nature of the model, the 
generating function G(t) of the statistics of the work 
is the product of those {Gk(t)}k corresponding to each 
mode with momentum k, i.e., G(t) = YikGk{t). The 
function G k (t) for a single harmonic mode is calculated 
in Eq. (All) of appendix A and is given by 



G k (t) 



1 



1-A 



-2iujf i t ' 



(15) 




FIG. 4. Probability density P(w) of the intensive work w = 
W/L d done on the one-dimensional [d = 1) bosonic model in 
Eq. (3) with m = 1, mo = 0.1, and of size L increasing from 
150 to 500 in steps of 50. For large L the prominent feature 
is a bell-shaped distribution which is completely described by 
its mean value, standard deviation and skewness (in order of 
decreasing importance) and which approaches a gaussian as 
L increases. The dashed vertical line indicates the asymptotic 
value of the mean for L — > 00. 



Accordingly, for L — »• +00, 



lnG(t) 



-ifb t + 



1 



ln(l 



In (1 - X 2 k e~ 2iuikt 



A 2 fe ) 



(16) 



where f b = lim^oo AE gs /L d = J k (uj k - u ok )/2. The 
first term on the r.h.s. can be omitted if, consistently 
with the conventions introduced above, we rename W as 
W + AE gs so that we measure W starting from AE gs . 
Note that the integrals involved in the second and third 
term on the r.h.s. of Eq. (16), with ui k = y/k 2 + m 2 , are 
separately convergent in the ultraviolet (UV) k — > 00 pro- 
vided that d < 4 because X k oc k~ 2 in the same limit. By 
comparing Eq. (16) with the decomposition in Eq. (8), 
one identifies the second term on the r.h.s. with — 2f s 
while the third one turns out to be the analytic con- 
tinuation of the critical Casimir contribution —fc{R) hi 
Eq. (7) to imaginary film thickness R — > it. We consider 
this issue in more detail in section. IV 

The cumulants of P(W), their scaling with the system 
size L, the fidelity J 71 / 2 , and the profile of the lowest edge 
singularity can be derived from Eq. (16). In particular, 
the cumulants n nW can be calculated from lnG(i) using 
the relations 



i n — lnG(t) 



(17) 



The resulting expressions for K n w involve an integration 
over the momentum k, which we have evaluated in the 
continuum limit, i.e., J, = f2 d / (27r) d J °° dk k d ~ x when- 
ever the corresponding integral with a large-momentum 
cutoff A (related to the inverse lattice spacing on a lat- 
tice) converges for A — > 00. Depending on the space 
dimensionality one finds for the mean work (i.e., the so- 
called irreversible work ): 



2 \2 



(mg — m ) 



L" 



k 4w fc(w/c +uj 0k ) 2 (2n) c 



-X 



[mg — m 2 + 2m 2 ln(m/m )]/4 
(mo — m) 2 (mo + 2m)7r/6 
[(m 2 — 3m 2 ,) (to 2 — to 2 )/4 — to 4 ln(??i/ Too) 
+ (to 2 , - to 2 ) 2 ln(2A/TO ))]7r/4 



in Id, 
in 2d, 

in id. 

(18) 



Notice that nw is UV-convergent in d = 1 and 2 and 
it is generically finite upon quenching to the critical 
point to = 0. On the other hand, it is finite for 
d = 2 and 3, when the initial state approaches a crit- 
ical one too = 0, while it diverges in the same limit 
for d = 1. Note that the mean work (W) can be ex- 
pressed as (^>o\H(m) — H(m )\^>o), i-e., from Eq. (3), 
(W) = (to 2 - m§)(#o|/ fc <£k</>-fc/2|# }, where this lat- 
ter expectation value depends only on too and therefore 
(W) depends linearly on to 2 : the additional dependence 
of nw on to displayed in Eq. (18) comes about because 
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W is measured here with respect to AE gs = £7(0) — Eq(0) 
(i.e., fiw is actually the mean irreversible work) and E(0) 
depends on m 2 in a non-linear fashion. Further below we 
will extend this argument to a thermal quantum quench 
in order to argue that (W) depends linearly on m 2 also 
in that case. 
The variance is 



~2 _ „ _ TO 



(m 2 



;(m,Q - to 2 ) 2 X 



(2tt)< 



7r/(8m ) 
- ln(m /A)7r/4 
-m 7r 2 /4 + 7rA/2 



in Id, 
in 2d, 
in 3d. 
(19) 



a^y is UV-convergent only in d = 1. Note that the rela- 



'w 

tive standard deviation is 



L -d/2 



(20) 



i.e., the relative fluctuations of W vanish for L 
Analogously, from Eq. (16) one finds 



k 3w = L a (m% - to 2 ) 2 



J 0k 



L« 



k 8a& " (27TY 

[{m/m ) 2 - 1 + 21n(2A/m )]/4 
(m — 

(3m 2 - 2m 2 )7r/4+ 

(to 2 - 2to 2 ,)(tt/2) ln(2A/m ) + ttA 2 /2 



(m 2 , 



to 2 ) 2 x 



in Id, 
in 2d, 

in 3d. 
(21) 



The skewness of P(W) is defined as 



7i w 



K 3W 



'W 



-d/2 



(22) 



i.e., it vanishes for L — > +oo. As heuristically expected 
on the basis of the central limit theorem, the expres- 
sions above show that the mean work scales proportion- 
ally to the system volume L d , while the relative vari- 
ance and skewness are proportional to the inverse square 
root of the volume, indicating that the probability distri- 
bution P(W) of the intensive work w (or, alternatively, 
of the normalized work Wj (W)) becomes asymptotically 
sharply peaked as L increases. In fact, the scaling of the 
distribution P(W) for large L can be easily derived: ac- 
cording to Eq. (1G), lnG(£) is proportional to L d , i.e., it 
is extensive and therefore each of the cumulants K nW is 
extensive too. It is more interesting to study the prob- 
ability distribution of the work per unit volume (inten- 
sive work) w = W/L d whose characteristic function is 
G(t) — G(t/L d ). As long as one is interested in small rel- 
ative fluctuations of the intensive work around its mean, 
i.e., in the behaviour of P(W) close the central peak (see 
Fig. 4), it is possible to expand the function \nG(t)/L d 
around t = 0: 



\nG{t)/L d = l-iat- 



\bt 2 + ±ict 3 +C(i 4 ) (23) 



where the constants a, b, c become independent of L for 
L large and are given by a = fiw/L d , b = / L d and 
c = K-3w/L d - Accordingly, the asymptotic form of G(t) 
is 



G(t) =exp 



-iat 



1 

21? 



(b-a 2 )t 2 + 



3lL 2d 



(c + 2a 3 - 3ab)t 3 + O 



(24) 



which is governed by the lowest-order cumulants. More 
specifically, the probability distribution of w for large L 
is gaussian with mean value a and variance (b — a 2 )/L d 
which vanishes as L — > +oo. Figure 2 shows typical 
distributions of the work P(W) for various values of L 
while Fig. 4 displays the asymptotic shape of P(w) as L 
grows. 

In order to characterize the behaviour of the probabil- 
ity P(W) it is useful to calculate the square fidelity F 
between the pre- and post-quench ground state. Recall- 
ing its definition in Eq. (12), T is given by the exponential 
of the second contribution in Eq. (16), i.e., 



^(toq, to) = exp 



Jk VV^o/c^fc 



(25) 



Note that T(mo,m) = F(m, mo). The integral is UV- 
convergent for d < 4, in which case it can be expressed 
in terms of elliptic integrals and — In J- takes a scaling 
form of to and too In particular, — (moL)~ d In T is the 
function of to/too plotted in Fig. 5 for various values of 
d and which, for m — (or, equivalently, for too — > oo) 
takes the values 



Id, 



lnJ"(TO O ,0) = - ( : { 7r(21n2 - l)/4 in 2d. 

7r(37r- 8)/9 



d [2(l-7r/4) 

' ) , , . .„ 

in 3d. 

(26) 

For later convenience, we report here explicitly the be- 
haviour of J(mo,m) in the two limiting cases for the 
value of the masses to and too: 

.F(to ,to)~^ ^ {mL ) d c, ( ^ ( 27 > 

I e y > ° d for too -C to, 

where Cd is the value for to/to = of the function 
plotted in Fig. 5, which can also be inferred from Eq. (26). 

As we have anticipated, some of the expressions in 
Eqs. (18), (19), (21), and (25) for the first three cumu- 
lants and the fidelity, respectively, or their generalized 
susceptibilities (i.e., their derivatives with respect to to 
or mo) diverge for m — > and/or mo —> 0, especially 
for d = 1. These infrared (IR) divergences indicate the 
emergence of a universal scaling behaviour of the cor- 
responding quantities in the neighbourhood of critical 
points. In particular for to — > 0, the mean work sus- 
ceptibility d 2 ^,w/9m 2 (x ln(m/m ) exhibits a logarith- 
mic IR-divergence in Id for m or mo — > (see Eq. (18)). 
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FIG. 5. Scaling plot of the square fidelity T (see Eqs. (12) 
and (25)). The dimensionless quantity — (moL)~ d In .7- is plot- 
ted as a function of m/mo < 1 in d = 1 (red), 2 (green) and 3 
(blue curve). Due to the symmetry of T under the exchange 
m o mo, the behaviour for m/mo > 1 can be readily inferred 
from the present plot. 



IV. THE EDGE SINGULARITY AND ITS 
RELATION TO THE CRITICAL CASIMIR 
EFFECT 

The behaviour of P(W) for W close to the threshold 
is determined by the large-i behaviour of G(t), which, in 
turn, is given by the third contribution in Eq. (16): 



) = m 



(28) 



In order to calculate this asymptotic behaviour, we first 
expand the logarithm as a Taylor series, recalling that 
|Afe| < 1 for all finite values of to, Too so that the expan- 
sion is always justified. Next we interchange the order 
of momentum integration and summation over n, since, 
recalling that A^ oc k~ 2 as k — > +oo, all intermediate 
expressions are UV-convergent at least as long as d < 4, 
and therefore 



A k e 



(29) 



Analogously, the fidelity J 71 / 2 develops a non-analytic 
behaviour upon approaching the critical point m = 0. 
Indicating by g the linear coupling of the most relevant 
operator which drives the system away from the critical 
point (</> 2 in the present case — see Eq. (2) — and there- 
fore g = m 2 and go = to 2 ,), one defines generalized fidelity 
susceptibilities as Xn (go, g) = —L^ d dg\nJ 7l / 2 (g ,g) 38 ' 39 . 
In particular, the fidelity susceptibility X2{g,g) is ex- 
pected to scale as x 2 {g,g) ^ \g ~ g c \ ud ~ 2 for g -t g c 
(see, e.g., Ref. 14, 38, and 39), where g c is the critical 
value of the coupling g and v is the bulk critical ex- 
ponent which controls the divergence of the correlation 
length £(g —> g c ) ~ \g — g c \~ v of the system. In the 
present case g c = and, as expected for a gaussian the- 
ory, v — 1/2. A direct calculation from Eq. (25) shows 
that X2(m- 2 ,™ 2 ) oc (m 2 ) d / 2 ~ 2 for to — > and d < 4, 
in agreement with the expected scaling. For m — > 0, 
instead, IR-divergences are even more common and ap- 
pear also at the level of the cumulants themselves, as one 
can easily see from Eqs. (18), (19), and (21). A physical 
interpretation for this is that the quadratic potential of 
the pre-quench hamiltonian is vanishingly small for the 
modes with low wavevectors k (as u;fe_ > o( m o = 0) — k), 
the wavefunction of which is therefore extended in space, 
such that an infinite amount of work is required in or- 
der to raise the potential to its post-quench form with 
<^fc->o("^) = fn- This fact has an important consequence 
also for the statistics of the large deviations of the work 
W from its mean (and typical) value, and leads to the 
occurrence of a condensation transition . 



If to ^ 0, the stationary phase approximation can be 
used in order to determine the behaviour for rat 3> 1 of 
each integral in the series with u>k — rn + k 2 /(2m), which 
therefore decays as 



\ In —InibJk t 
A k e 



\2n 
A 



' TO \ d / 2 



(30) 



with Ao = (too — to) /(too + to). The series in Eq. (29) 
can then be summed up and it gives 



<p(t) 



1 / m \d/2 _ 

2 \M) L1 W(V 



(31) 



where lA s (x) is the polylogarithm function of order s (see, 
e.g., §25.12 in Ref. 40). 

Recalling from Eqs. (9), (16), (12), the discussion after 
Eq. (16), and Eq. (28), that the work distribution P(W) 
is given by 



P{W) =J r (m ,m) 



dt 
2?r 



jWt-L d v {t) 



(32) 



one can see (by formally expanding the exponential in ip) 
that each term in the Fourier series of the (7r/TO.)-periodic 
function Li^^^oe^ 2 "™*) results in an edge singularity 
in P(W) located at W = 2nm, with n = 1,2, There- 
fore the lowest threshold of P(W) is at W = 2m and the 
associated edge singularity is determined by the lowest 
term of such a Fourier series, and — as discussed further 
below — it takes the form oc d{W - 2m) (W - 2TO) d / 2 " 1 , 
where ^{x) is the Heaviside step function. This heuristic 
argument is supported by the following analysis, which 
allows a calculation of the expression of the probabil- 
ity P(W) for W < 4to beyond the stationary phase ap- 
proximation. First of all we observe that if the (inverse) 
Fourier transform f{ui) of a function f(t) vanishes iden- 
tically for ui < ojq, then the Fourier transform of e^** 1 is 
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given by 6(v) + f(w) for u < 2uj q . In order to see this, 
one can formally expand e^*) = 1 + f{t) + f 2 (t)/2 + ... 
and note that, by virtue of the convolution theorem, the 
n-th term of this expansion contributes to the Fourier 
transform of e-^'J only for u > tiljq. Accordingly, for the 
specific case we are interested in, P(W) for < W < 4m 
is proportional to the inverse Fourier transform <p(W) of 
tp{t), which, in turn, is given only by the inverse Fourier 
transform of the first term of the series in Eq. (29), since 
all the other terms contribute to (p{W) only for W > 4m. 
Accordingly, P(W) for < W < 4m, i.e., the profile P es 
of the lowest edge singularity, is given by 



P cs (W)= 1 -L d F( mo ,m) J ^j Wt j 



A 2 e -2^ fc t 



l -L d F{m ,m) [ \ 2 k 6(W - 2co k ) 



2 

(mL) d 



J"(mo,m) x ( jr— ; — ) > (33) 
in \ 2m m ) 



with 



x x 



\J x 1 + p 2 — 1 — X 
\/ x 2 + p 2 — 1 + x 



(34) 



Note that the first term of the expansion of the exponen- 
tial contributes to P(W) with J r 5(W), as expected. For 
x ~ 1 + , i.e., upon approaching the threshold W — >• 2m + , 
Eq. (34) becomes 



X(x -> l + ;p) 



4(2^; 



- d \l[2{x-l)f 2 -\ (35) 



[where A = (p — l)/(p + 1)] which renders the algebraic 
behaviour of P(W — > 2m) that we anticipated above and 
that can be alternatively inferred by taking the inverse 
Fourier transform of the function tp(t) calculated within 
the stationary-phase approximation mt ^> 1 reported in 
Eqs. (30) and (31). In deriving these results, we assumed 
that the system size is large but still finite, otherwise 
the expansion of the exponential e~ L v ^> is no longer 
justified and the problem can be studied only in terms 
of large deviations of the intensive work w — W/L d , as 
discussed in Rcf. 14. 

The behaviour of P(W) at the threshold is generically 
non-analytic: when the threshold is approached from 
above, P(W) diverges for d < 2, attains a finite value 
for d = 2, while it vanishes for d > 2. Figure 3 shows 
(for a choice of mo and m) the profile of the lowest edge 
singularity in d = 1, 2, and 3 as obtained numerically 
by a direct evaluation of the integral in Eq. (32), with ip 
and T given in Eq. (28) and (25), respectively. This nu- 
merical result (solid line) is compared with the analytic 
expression P cs (dashed line) provided by Eqs. (33) and 
(34). The agreement between these curves is excellent 
for W < 4m, as expected. On the other hand, one can 



easily verify that the asymptotic expression for W —> 2m 
in Eq. (35) — determined by the behaviour of ip(t) for 
mt 3> 1 — does not provide an accurate approximation 
of P cs as soon as one moves away from the threshold, as 
it is clear by comparing Eqs. (34) and (35), which dis- 
play a qualitatively different dependence on x and differ 
already at the linear order in x — 1. In order to obtain 
a better approximation of the actual dependence on x of 
the scaling function x( x ', p) it is necessary to account for 
the behaviour of ip(t) also for mt ~ 1. Further below we 
discuss how to extend the analysis of the scaling limit 
to generic values of mt, highlighting the emergence of a 
scaling function which is expected to fully describe it. In 
particular, the connection with the critical Casimir effect 
that we establish below for the present non-interacting 
model is a special case of the more general correspon- 
dence discussed in Ref. 12 — and exemplified therein for 
quenches in the quantum Ising chain in a transverse field 
- which allows one to determine the behaviour of the 
edge singularity in interacting cases. Notice, however, 
that one generically expect that corrections to the scal- 
ing limit — which depend, e.g., on the lattice spacing 
a and on additional microscopic parameters — become 
increasingly relevant for values of W and 2m which are 
large on the microscopic scale. In this case, the behaviour 
of the system and therefore the statistic of the work is no 
longer solely determined by the long-wavelength modes 
in terms of which universality typically emerges. 

The discussion in the previous paragraph assumes that 
G(t) for large L is given by Eq. (16), which is the 
case only if the sums L~ d J^k over the wave vectors 
k n oc n/L, n 6 Z d allowed in a system of finite extent 
L (and suitable boundary conditions) can be replaced 
by the integrals J, indicated on the r.h.s. of Eq. (16). 
In particular, the discreteness of the set of allowed mo- 
menta k n implies that the lowest excitation energies 2u>k n 
of each single harmonic oscillator form a discrete set and 
that, accordingly, W takes values in this set (see, e.g., 
the second line of Eq. (33) with the integral replaced by 
the sum). For a given gap m of the post-quench hamil- 
tonian, the difference between the possible subsequent 
values of W (related to the inverse density of states) de- 
creases quickly as W increases above the threshold m. 
However, the first allowed value above the threshold 2m 
is given by 2wfe 1 with kf (x 1/L 2 . In order for the spac- 
ing AW = 2uk 1 — 2m to be negligible compared to the 
scale of variation of P(W), set by the distance 2m be- 
tween two subsequent thresholds (see the scaling form in 
Eq. (33)), one has to require AW <C 2m and therefore 
L 3> m _1 , i.e., the size L of the system has to be always 
larger than the inverse gap of the post-quench hamilto- 
nian or, equivalently, to the correlation length £ within 
the system. This condition is never fulfilled at the criti- 
cal point m = 0, which is characterized by the fact that 
P(W) varies on a scale set by ~ 1/L and therefore the 
discrete set of allowed values of W cannot be replaced 
as a continuum one close to W ~ 0. In addition, the 
consecutive thresholds identified above merge into a sin- 
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gle one at W = and the argument which allows one 
to relate P(W) to the inverse Fourier transform of tp(t) 
does no longer apply because all the subsequent terms of 
the expansion of the exponential, which are higher-order 
powers of ip(t), contribute to P(W). For these reasons 
we consider below only the case in which the inverse gap 
mT 1 is large on the microscopic length scale set, e.g., by 
the lattice spacing a, such that the continuum limit is 
justified, but it is still significantly smaller than the size 
L of the system. 

Following the general discussion done after Eqs. (7) 
and (8) about the relation between the generating func- 
tion of the work statistics at a quantum quench occurring 
close to a critical point and the critical Casimir effect, 
here we analyze this connection in more detail. In partic- 
ular, we would like to compare the form taken by ip(t) (see 
Eqs. (28) and (16)) in the scaling limit with the universal 
Casimir contribution to the free energy of a classical fluc- 
tuating gaussian (quadratic) field (equivalent to a system 
of quantum non-interacting bosons) confined in a film of 
finite thickness and with suitable boundary conditions. 
It is well-known that the universal behaviour of this sys- 
tem upon approaching the critical point can be described 
in terms of a small number of effective boundary condi- 
tions (or, more precisely, in terms of the so-called surface 
universality classes 12 ' 41 ), depending on some gross fea- 
tures of the quench under consideration (e.g., breaking 
of some symmetry etc.). In the context of quenches in 
quantum systems ' it has been argued and then shown 
explicitly for non-interacting bosons that the proper ef- 
fective boundary condition for the field is of Dirichlet 
type (see also section III.B in Ref. 42 for further discus- 
sion), at least for sufficiently "deep" quenches, i.e., with 
an initial mass m which is large compared to all the 
other scales at play in the system (i.e., m and possible 
microscopic scales). As we will see below, our comparison 
also support this conclusion. 

We first show that <p(t) in Eq. (28) acquires a scaling 
behaviour for t and m large enough on a microscopic scale 
(but, contrary to the analysis leading to Eqs. (30), (31), 
and (35), with mt assuming generic values). In order 
for the integral in Eq. (28) to be well-defined also for 
Xk — > 1, we will assume the prescription t — > £ — z0 + , 
which corresponds to introducing an ultraviolet cutoff in 
the momentum integrals. Then, in this scaling limit, each 
integral in Eq. (29) is determined by the low-fc behaviour 
(k <C a -1 where a the lattice spacing) of its integrand 



\2n — 2niWfct 
A k e 



A,-," ,"'!, 2 "'* / d^ (w 2 - m 2 ) d/2 e- 2niut , 



° (2TT) d d 



(36) 



where an integration by parts has been applied, in which 
the contribution for w& — > oo vanishes due to the pre- 
scription t = t — i0 + . In Eq. (36), Xk has been approxi- 
mated by Ao on account of the fact that a direct inspec- 
tion of its expression (see Eq. (14)) shows that Xk as a 
function of ujk varies on a scale set by mo, whereas the 



exponential factor varies on a scale ui ~ Accord- 
ingly, for \t\ 3> Mq 1 the variation of the former can be 
neglected. Note that in principle one could account for 
the complete fc-dependence of Afe in Eq. (36), as it is done 
further below when discussing the case mo = 0; however, 
for the purpose of highlighting the connection with the 
critical Casimir effect it is more convenient to approxi- 
mate first Afe ~ Ao- Now summing the series Eq. (29) 
and using Eq. (36), we obtain 



(p(t) = (it)~ d Q(imt;X Q ), 



(37) 



where 



- 1' 



A plot Q(x; A) as a function of x is provided in Fig. 6 for 
various values of A and d = 1. The scaling form (37) for 
tp{t) can be directly compared with the known results for 
the critical Casimir effect in a film. Note that the model 
(2) we are presently interested in involves a quadratic 
hamiltonian of the real scalar field <f> associated with the 
relative phase of the quasi-condensates and therefore its 
classical correspondent is a real scalar field, correspond- 
ing to the case N — 1 of the iV-component scalar field 
with O(N) symmetry discussed in Ref. 43. The univer- 
sal Casimir contribution fc to the free energy per unit 
area S = L D ~ l and per temperature /3" 1 of such a fluc- 
tuating classical gaussian (quadratic) field in D spatial 
dimensions confined in a film of thickness R with Dirich- 
let boundary conditions (corresponding to the so-called 
ordinary surface universality class ,44 ) is given by 1 



fc(R,(;D) = R-^@ o (R/0, 



(39) 



where £ is the correlation length of the fluctuations of 
this classical gaussian field in the bulk (i.e., for R —> oo) 
and 



,.D 



e o (a0 = - 



( 47r )(D-l)/2 r (£±1) j 1 



1 



(40) 

is the associated universal scaling function; see, e.g., 
Eq. (6.6) of Ref. 43. By comparing Eqs. (37) and (38) 
with Eqs. (39) and (40) one concludes that <p(t) can be 
cast in the form 



<p(t) = f c (R = it + + ,S 



D 



I) 



(41) 



for Ao — > ±1, i.e., either mo 3> m or mo <C m. However, 
we shall see below that the latter case actually requires 
accounting for corrections to Afe which have been ne- 
glected in deriving Eq. (36). Equation (41) demonstrates, 
for this specific case of the bosonic field, the connection 1 2 
between the generating function G(t) and the critical 
Casimir effect free energy fc(R) of the corresponding 
classical model, which we discussed in section II. As an- 
ticipated above, the scaling function Q(x, A) which con- 
trols the behaviour of fit) (see Eqs. (37) and (38)) coin- 
cides with the one 0o of the critical Casimir force in the 
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presence of Dirichlet boundary conditions (see Eq. (40)) 
only in the limit mo > mof a deep quench. In this case, 
the overall amplitude of the edge singularity, which is set 
by T according to Eq. (33), is exponentially suppressed 
as a function of (moL) d Cd (see Eq. (27)), where Cd gener- 
ically takes rather small values. Accordingly, in order for 
this feature to be observable one has to be able to de- 
tect somewhat rare events. Furthermore, based on the 
discussion done after Eq. (31), the behaviour of P(W) 
for 2m < W < 4m is generically determined by the in- 
verse Fourier transform of —L d ip(t) and therefore of the 
associated universal scaling function 8 (see Eq. (37)). In 
particular, from Eqs. (37) and (38) one finds that P(W) 
takes the scaling form of Eq. (33) with a scaling function 
X 1 given by 



4(2tt 



Ana; 



(* 2 - 1) 



d/2-l 



(42) 



for 1 < x < 2 and An = (p — l)/(p + 1). Compared 
to the scaling function in Eq. (35), which accounts only 
for the behaviour of (p(t) for mi 3> 1 and it is therefore 
accurate only for x — > 1 + , yf provides a better approx- 
imation of the actual scaling function in Eq. (34) and 
for p ^ 1 reproduces it in the full range 1 < x < 2. 
Note that, according to this relation between P(W) and 
Q(x; A), the leading behaviour of P(W) upon approach- 
ing the threshold, i.e., for W — > 2m + depends on the 
behaviour of the scaling function at large (imaginary) 
values of the scaling variable mi, which is connected to 
the asymptotic behaviour of /c(-R) a t large (real) val- 
ues of R. In fact, rather generally, the relation between 
P(W) and G(t) in Eq. (9) implies that G(t) as a func- 
tion of the complex variable t = tji + itj G C, with 
</.7V, € K, cannot have poles in the lower complex half 
plane ti < because P(W < 0) = (here we assume 
that W is referred to the threshold value AE gs ). In ad- 
dition, G(t) has the symmetry G*(t) = G(-t*). For 
the specific model we are presently interested in, G(t) in 
a finite volume (e.g., assuming periodic boundary con- 
ditions) is the product over the allowed wave-vectors k 
of the generating function Gk of each single mode, i.e., 
G (t) = lL G >c(0! with G k {t) given in Eq. (15). Each 
of these factors is characterized by a branch cut when- 
ever the argument of the square root (or, equivalently, of 
the logarithm in lnG(t) oc ip(t) in Eq. (16) before taking 
the limit L — > oo which turns into L d f.) becomes 
real and negative, i.e., for tj > (In A^ 2 ) / (2wfc) > and 

tu = t-n = mr/(2cuk), with n 6 Z. (In passing, we 
mention that the analogous singularities for the quan- 
tum Ising chain in a transverse field are represented by 
simple poles .) For a given mode fc, these cuts are par- 
allel to the imaginary axis, include it, and are displaced 
along the real axis. As expected, however, they are all 
located in the upper half plane ij > 0. Accordingly, G(t) 
and (p(t) are analytic in the lower half plane t j < and 
therefore the behaviour of ip(t) at large t along the real 
axis (with a vanishingly small negative real part — i0 + ) 
can be obtained by analytic continuation of the corre- 



sponding behaviour of the function along the negative 
imaginary axis —iR with R > 0, where it coincides with 
the critical Casimir force fc acting within a film of thick- 
ness R, as anticipated in Eq. (41). 




FIG. 6. Scaling functions Q(x; X) (solid blue curves) for 
A = 1, 0.9, 0.8, 0.6, 0.2 and @ c (x) (dashed red curve) for the 
non-interacting bosonic model in one spatial dimension (see 
Eqs. (38) and (47) respectively). In contrast to Q(x; A), O c {x) 
vanishes for x — 0. 

The results discussed above are valid in the scaling 
limit and for \t\ 3> , which allowed the approxima- 
tion Xk — Ao in Eq. (36). The leading correction to this 
approximation takes the form of a rcnormalization of t 
and Ao- In fact, the first-order correction to A& as a func- 
tion of uik can be expressed in exponential form 

A fe = Aoe- 2 ^"™)/™ [1 + 0(( Wfc - mf/ml)] (43) 

and therefore it translates into an imaginary shift of t 
and a rescaling of Ao in Eq. (37): 

t^ i cff = t-2i/m and A ^ A cff = A e 2m/m °, (44) 

i.e., to the introduction of the effective parameters t e s 
and A e ff. In particular, the emergence of i e ff confirms 
explicitly the expectation 1 that a large but finite value 
of mo amounts at an extrapolation length 



(45) 



at each of the two surfaces of the film in the bound- 
ary formulation of the quantum quench problem; £ ext ac- 
counts for the fact that the fixed-point boundary condi- 
tion (in this case of Dirichlet type for mo large enough) 
is not effectively imposed at the boundaries of the film 
but on displaced effective boundaries. Analogously, the 
parameter A e ff accounts for the effective "strength" of the 
boundary condition and it approaches 1 for mo 3> m, i.e, 
p = mo/m 1. The actual effect of mo taking a finite 
value > m (i.e., p ~ 1) on the behaviour of the probabil- 
ity distribution P CS (W) for 2m < W < 4m is described 
by Eq. (33) with the scaling function x in Eq. (34) , which 
was derived by considering the full fc-dependence of A^ in 
Eq. (29), while it is not entirely captured by the scaling 
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function X ' m Eq. (42) which corresponds to £ cx t = 0. 
However, part of this effect can be accounted for by tak- 
ing the inverse Fourier transform of Eq. (37) with the ef- 
fective parameters in Eq. (44). Because of the translation 
t i— > £ e ff = t — 2i£ cxt , this transform acquires an overall 
multiplicative factor e~ 2Wimt compared to the case with 
^ext = while A i-> A e ff = A e 2 / p according to Eq. (44). 
The resulting P(W) takes the same scaling form as in 
Eq. (33) with a scaling function \" given by 



= e -4(x-l)/pV 



(46) 



With this additional exponential factor, the difference be- 
tween x" and the actual scaling function x turns out to be 
of order 0([(x — l)/p] 2 ), with a significant improvement 
compared to the scaling function x' ■ I n summary, while 
the knowledge of the asymptotic behaviour for large sepa- 
rations of the critical Casimir force in the classical system 
is sufficient to predict the leading behaviour of the prob- 
ability P(W) very close to its lower threshold, the knowl- 
edge of the full scaling function 8 provides the necessary 
information to predict P(W) with a certain accuracy up 
to the second threshold, while additional improvements 
can be obtained by considering the emergence of an ex- 
trapolation length and of an effective surface coupling 
within the corresponding classical system. 

Upon decreasing too, the extrapolation length £ cxt — 
TOq" 1 which controls the corrections mentioned above 
grows, the influence of the corrections increases and even- 
tually one can no longer neglect the full fc-dependence of 
Afe, which is particularly strong for too = 0. In this case, 
A/c(wfe)| TOo= o has to be accounted for explicitly in Eq. (36) 
where it amounts to an additional multiplicative factor 
that is integrated by parts in the next step. As a result, 
one finds the same scaling behaviour as in Eq. (37), but 
with a novel scaling function given by 



e c (a:) 



d(27r)° 



d d 
X 



As 



(s 2 - \) d ' 2 



(Vs 2 - 1 + s) 4 e 2xs 



1 



(47) 



Differently from the function in Eq. (38), the present 
scaling function vanishes for x — > 0, i.e., if the post- 
quench hamiltonian approaches the critical point, with 



Q c (x -> 0) = 



d r(i-f) 

C 2d(47r) d / 2 

4 r(d-4) 
C 2 d (i-x) d / 2 T(d/2) 



for d < 4, 



for d > 4. 



(48) 

This behaviour is actually expected since for to — > — 
Too the quench becomes vanishingly small and therefore 



no work is done on the system and the edge singular- 
ity disappears. Figure 6 shows plots of the universal 
scaling functions Q(x;X) and c (x) for d = 1. Note, 
however, that the functional form of the behaviour of C 
for x 3> 1 is the same as the one of 9(x;A ^ 0) and 
therefore the leading algebraic behaviour of the edge sin- 
gularity as W — > 2m in the case to = is the same 
as for too ^ 0, up to an overall amplitude, in agree- 
ment with Eq. (35). The inverse Fourier transform of 
the critical Casimir contribution (p(t) = (it)~ d Q c (imt) 
with 6 C given in Eq. (47) can be calculated and the 
probability distribution P(W) for 2m < W < Am ac- 
quires the scaling form in Eq. (33) (with mo = 0), with 
the scaling function X replaced by X c(%) = x{ x i 0). This 
scaling function for x —¥ 1 + has the same behaviour as 
X, x' an d x" U P to an overall multiplicative factor, be- 
cause the asymptotic behaviours of Q(x; Ao) and Q c (x) 
for x 3> 1 are the same. However, depending on the ac- 
tual spatial dimensionality d of the system, this common 
behaviour might emerge only very close to the threshold 
x = 1 because of the presenc e of rathe r strong corrections 

Xc(x)/X"(x,p) = Ao 2 [l - - 1) + °( x - !)] which 

mirror those in the asymptotic behaviour of <d c (x 1) 
compared to Q(x 3> 1, Ao) clearly visible in Fig. 6. More 
pronounced differences emerge upon moving further away 
from the threshold, which are connected to the markedly 
different behaviours of the corresponding scaling func- 
tions for x < 1. Note that in the present case mo = 0, the 
overall amplitude J"(mo = 0, m) of P{W) (see Eq. (33)) 
is exponentially decreasing as a function of (mL) d Cd, as 
indicated by Eq. (27), where Cd generically takes rather 
small values and mL has been assumed to be large com- 
pared to 1 in order for the present analysis to be valid. 
Accordingly, in order for the edge singularity P(W) to 
be observable one has to be able to detect somewhat rare 
events. 



V. THERMAL QUANTUM QUENCH 

In this section we extend our investigation of the statis- 
tics of the work to the case in which the initial state 
is not the ground state of the pre-quench hamiltonian, 
but a thermal state at temperature /3 _1 . The analy- 
sis which follows is an extension of the one presented 
in Ref. for a generic quench of a single oscillator to 
the case of a collection of non-interacting harmonic oscil- 
lators (see Eq. (3)) subject to an instantaneous quench. 
Due to the absence of interaction, the generating function 
Gth(i) of the thermal P(W) is still given by the product 
Gth(i) = life Gth,k{t) of the generating functions for each 
single mode k, which is calculated in Eq. (A24) of the 
appendix A. The result is 
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if b t+ / [ln(l 

1 



) + -ln(l-A 2 )]- 



ln 



/j _ g-2(/3— it)« fc^2\ _ 2^ — ^2\ e -(^— it)«o*— iw*t _ ^2 _ Q -2(J3-it)uj ok ^ e ~2iuj k t 



(49) 



Consistently with the conventions introduced after 
Eq. (16), the term —ifyt on the r.h.s. above can be omit- 
ted if one measures W starting from the difference AE gs 
between the energies of the ground states of the post- and 
pre-quench hamiltonians by renaming W as W + AE gs . 

We repeat the calculation of the cumulants of the dis- 
tribution P(W) for the thermal case and derive the lead- 
ing corrections they induce to the ordinary quench results 
in the limit of small temperatures (/3 — » oo). In particu- 
lar, for the mean work /i' w = iG' th {0) one finds 



(J-'w = Mw - L d (ml 



(50) 



ik 2w ofc (l 

where = iG'(0) is the work done during the corre- 
sponding ordinary quench = 0), given in Eq. (18). 
As expected, the difference (j,' w — fiw vanishes at low 
temperatures j3mo 3> 1 as 



-(m L) c 



m 



-/3m 



2m (2nPm ) d / 2 ' 
Similarly, for the variance a'w one finds 



° W 



+ L d (m 2 -m 2 ) 2 



k 2u: 2 k (l~e~^ok)2- 



(51) 



(52) 



where crj^ is the variance in the corresponding ordinary 
quench, given in Eq. (19); The lowest-order correction at 
low temperatures is 

s -/3m 



w 



>w 



(m L) 



2m 2 



(2irf3m ) d / 2 ' 



(53) 



These expressions show that thermal corrections are 
gcnerically of order ~ e _ ^ m ° . 

As in the case of the ordinary quench, it is interesting 
to compare the work done during the thermal quench 
with the work W rcv required in order to change m from 
its pre- to its post-quench value in a reversible way, i.e., 
while keeping equilibrium with the bath at temperature 
/3 _1 . According to thermodynamics, the latter is given 
by the change 

AF = F (g) - Fp(g ) (54) 

of the free energy Fp(g) of the system, where we indicate 
generically by g the parameter of the quench. In the case 
we are presently interested in we identify g with m 2 and 
go = m 2 ,, while Fp(g) can be expressed, for large L, as 

Fp{g) = - r 1 lnTr{e-^)} 

(55) 



p- 1 ^ J In Z(w fc (m = Vfl)). 



where Z(oj) = J2™=i e 



-/3bj(n+l/2) 



[2sinh(/3u;/2)]~ 



is 



the partition function of a single harmonic oscillator with 
characteristic frequency u. The reversible work is there- 
fore given by 



W„ 



sinh[/3w fc /2] 



AF = /T 1 ^ / In 

l k smh[/3cjofe/2J 



(56) 



(where uj k = uj k (m = y/g) and u} k = w fc (m = y/go)) and 
it has to be compared with (W) — (J,' w , which renders 
AE gs reported after Eq. (16) for j3 — > oo. Taking into 
account Eqs. (50), (18), and the shift of W by AE gs (such 
that (W) = n' w + AE gs ), one eventually finds 



(W) = 



_9o Ld 



1 



1 



k uj ak tanh(/3w fc/2)' 



(57) 



Interestingly enough, these expressions for (W) and W rev 
are similar to those for the thermal quench of the Ising 
chain discussed in Ref. 17. In addition, Eq. (57) shows 
that the dependence of (W) on g = m 2 is linear, as in 
the case of the ordinary quench (see section III). In fact, 
according to the argument outlined after Eq. (18), 



(W) 



Tr{e-P Ho (H-H )} 
Tr{e-0*o} 

. ^ Tr{e-e H <>dHo/dgo} 
{9 ~ 90 > Tr{e~^o } 

(g -9o)Fp (9o), 



(58) 



and therefore the non-equilibrium work can actually be 
expressed in terms of the equilibrium free energy F^{go). 
With the proper identifications, the same relation holds 
for the case of the Ising model studied in Ref. 17 (where 
the strength of the transverse field plays the same role as 
g = m 2 here). Accordingly, the irreversible work W m 
can be expressed as 

W m = (W) - AF 



(9 - 9o)F' p {go) + Fp(g ) - F p (g) > 0, 



(59) 



where the last inequality — which is actually expected 
on the basis of thermodynamics — follows from the fact 
that F'J(g) < 0, i.e., Fp(g) is a concave function of g, as 
it can be verified from Eq. (55). 

The irreversible work W nl is responsible for an irre- 
versible increase of entropy ASi rl (g, go) = (3W\„ > 0. 
For a small quench with g = go + 8g, Eq. (59) yields 

AS iTI (g +Sg,g ) = -(Sg) 2 f3F^(go)/2+0(f3(Sg) 3 ), (60) 
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in which the only contribution on the r.h.s. comes from 
W lev — AF because (W) depends linearly on 5g. For a 
fixed (small) value of 5g, the irreversible entropy produc- 
tion as a function of the pre-quench value go (or, alterna- 
tively, of the post-quench one g) has been shown to dis- 
play a pronounced maximum upon approaching the crit- 
ical point (located at g = g^ ) in a thermal quench of the 
quantum Ising chain ' . This behavior has been traced 
back to the fact that, as the energy gap closes at critical- 
ity, most of the work done during the quench is dissipated 
in exciting the system and this dissipation is responsible 
for the emergence of an intrinsic irreversibility 1 ' quan- 
tified by AS'in-. Rather generally, we can infer from 
Eq. (60) that the emergence of such a maximum in A£i rr 
is however related to the behavior of F'J(g ) for g —> g^ 
and therefore it is de facto determined by the equilibrium 
free energy Fp(go) of the system at finite temperature /3 
upon approaching the (quantum) critical point go — > <?q . 
Both in the present case of the bosonic hamiltonian in 
Eq. (3) and for the Ising chain discussed in Ref. 17, the 
coupling constant g (and therefore go) couples the most 
relevant operator which drives the system away from the 
critical point and therefore g plays the role that the de- 
viation of the temperature T from its critical value T c 
has for a classical phase transition at finite tempera- 
ture. Accordingly, the possible non-analytic behavior in 
Fp(go) is characterized by the critical exponent a of the 
specific heat C oc —Fp(g ), i.e., C ~ \go — go\~ a an d 

therefore ASi„ ~ P(Sg) 2 \go — 5ol~ Q at l east as l° n g as 
\Sg\ <C \go — <7g|, such that higher-orders in the expansion 
in Eq. (60) are actually negligible. For the Ising chain 
at zero temperature = — which displays the same 
critical behavior as the two-dimensional classical Ising 
model at finite temperature — a = and the singularity 
is actually logarithmic. However, at finite temperature, 
there is no phase transition in the model and therefore 
the pronounced maximum which might be observed in 
A Si rl as /3 _1 — > is weakened, as it is clearly visible in 
Fig. 1 of Ref. 17. 

For the system with hamiltonian (3) we are presently 
interested in, we recall that the critical point is located 
at <?q = and go as well as g can take only positive val- 
ues in order for the hamiltonian to have a ground state. 
The corresponding free energy Fp(g) in Eq. (55) and in d 
spatial dimensions can also be expressed (after a suitable 
regularization of the UV behavior) in terms of the parti- 
tion function Z c \ cx / T)<\> exp{— S [(/>]} of a classical (real) 

field 4>(x, t) with gaussian action S = ^ dr J d d x C, 
where the (euclidean) lagrangian density £ is given by 
Eq. (2) after the replacement II i— )• d T (f>(x, r) and m 2 h->- g. 
The additional "spatial" dimension r is of finite extent 
/3 and the field 4> satisfies periodic boundary conditions 
4>{x,t = 0) = 4>{x,t = j3) along the r-direction, i.e., at 
the boundaries of this effective film of thickness (3. Note 
that this construction is analogous to the one discussed 
after Eq. (6) and therefore also in this case the free energy 



Fp(g) for large L decomposes as in Eq. (7): 

PF p {g) = L d [Pe b + f^ r \p)l (61) 

where £&(<?) = J k uJk/2 is the energy density of the ground 

state, whereas f^ er ^ (which depends, inter alia, on g) 
represents the finite-size contribution which vanishes for 
j3 — » oo and is responsible, close to the critical point, 
for the emergence of the universal critical Casimir ef- 
fect. Compared with Eq. (7) the "surface" contribution 
corresponding to 2f s is absent here because the periodic 
boundary conditions do not break translational invari- 
ance across the "film" "'. Upon approaching the critical 
point go — > <?o onc expects f(j er \l3) to take the scaling 
form in Eq. (39) where D = d + 1, £ = mr 1 = ,g~ 1//2 , 
R is replaced by f3 and the scaling function Qo( x ) ( see 
Eq. (40)) by the analogous function Q per (x) appropri- 
ate to periodic boundary conditions. Q per {%) is given by 
(see, e.g., Eq. (6.8) in Ref. 43) 

x D f°° {s 2 -!) 1 ^ 1 

&per{x) = ~(4 7T yD~i)/2 r{ D±l ) J 1 ds e --l ' 

(62) 

as one can verify by analyzing Eq. (55) after the subtrac- 
tion of the bulk contribution L d f3eb (see below Eq. (61)). 

Note that the specific heat C oc —F'^{g) can be conse- 
quently decomposed as in Eq. (61). In particular, upon 
decreasing the temperature, the effective thickness j3 of 
the film diverges and only the bulk term contributes to 
F^ (g) . Being the system effectively d + 1-dimensional, 
one expects C (and therefore F^(g) = L d e'^(g)) to be 
characterized by a critical exponent a = ad+i, where 
oljj = 2 — D/2 is the specific- heat critical exponent of a 
gaussian (scalar and real) field in D space dimensions. 
Accordingly, the expansion of £b(g) for small g contains 
a leading non-analytic contribution ~ g2-a d+1 j n ac j(jj_ 
tional to regular, analytic terms. At a finite temperature 
/3 _1 7^ 0, instead, the d + 1-dimensional system has a 
finite thickness: when the correlation length £ = g^ 1 ^ 2 
of the field fluctuations is much larger than the thickness 
f3 of the film — which occurs sufficiently close to the 
critical point — the system displays the critical behav- 
ior corresponding to d spatial dimensions, with a = ay 
and therefore one expects F'^(g) ~ g~ ad ■ In fact, from 
a suitable expansion for small argument of O per (x) in 
Eq. (62) one can verify that the resulting expansion of 

fc (P) in 9 f° r P 2 g *C 1 contains both /3-dependent 
analytic contributions (3~ d [ao + ai/3 2 <? + a2((3 2 g) 2 + . . .] 
and a leading /3-independent non-analytic contribution 
cx g d / 2 [l + 0(j3y/g)]. When inserted into Eq. (61), the 
latter term renders the expected singularity ~ g~ ad of 
F'I{g) as j — > 0, which indeed prevails on the non- 
analytic behavior of e' b ' discussed above. The presence 
of these non-analytic contributions follows from general 
properties of the scaling function &o,per(x) of critical 
Casimir forces and can be also verified via a direct cal- 
culation of F'J(g) based on Eq. (55). 
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FIG. 7. Probability density P(W) of the work for a thermal 
quantum quench in Id with m = 1, mo = 0.1, /3 _1 = 0.1, 
and increasing system size L — 50, 75, 100, and 125. Ex- 
cept for the non-zero temperature, all other parameters of 
the quench are the same as in Fig. 2, with which the present 
figure should be compared. Note that, in contrast to Fig. 2, 
there are peaks also for W < 2m which actually extend into 
the region W < (not shown here because these additional 
peaks are not visible on the vertical scale of the present plot 
due to their tiny height). All these peaks are of purely ther- 
mal origin and correspond to transitions originating from the 
thermal excitations of the initial state. 



Accordingly, the behavior that the irreversible entropy 
production AS lrr (go + Sg,go) displays for small Sg upon 
varying g and therefore the singularity developed upon 
approaching its critical point <?o = ffo = 0> is completely 
controlled by equilibrium critical exponents and a dimen- 
sional crossover 30 '' 51 — resulting from the competition 
between the bulk and the finite-size term in Eq. (61) - 
is observed by lowering the temperature f3~ , which is 
accompanied by an increase of A5i rr oc /3. 

Figure 7 shows P(W) for a specific thermal quench 
of the one-dimensional model in Eq. (2) and for vari- 
ous values of L. Differently from the zero-temperature 
case reported in Fig. 2 for the same values of the quench 
parameters but with /3 _1 = 0, the lower threshold at 
W = 2m has disappeared along with the regular struc- 
ture of an equally spaced sequence of peaks. In addi- 
tion, W can now assume arbitrarily large negative values 
(though with an exponentially small probability) due to 
the fact that the system can be in an arbitrarily highly 
excited state of the pre-quench hamiltonian, which can 
provide (W < 0), rather than absorb (W > 0), the en- 
ergy required for the formation of post-quench excita- 
tions. Note that these features of Fig. 7 are actually 
generic and their occurrence does not depend on the spe- 
cific choice of the parameters of the quench as long as 

According to the general discussion done in Ref. 46, we 
expect the distribution function P(W) generated by G t h 
in Eq. (49) to satisfy the fluctuation relations which char- 
acterize the statistics of the work done on isolated quan- 
tum systems. Indeed, in the present case of a quench, 



it is rather straightforward to show that the Crooks and 
Jarzynski fluctuation relations ' ' are satisfied, following 
essentially the analysis of Ref. 19. G(t) in Eq. (11) can 
be cast in the form 

ry r r p-iHt iH (t+i/3)l 

G Ho ^h((3; t) = (e- iW <) = ir|C 7 ° ] - (63) 

where Zh ((3) — Ti{e~^ H °} is the partition function of 
the system with hamiltonian Hq at inverse temperature 
(3. In this expression we have clearly indicated with a 
subscript that G refers to the quench from Ho to H . The 
cyclic property of the trace, implies 

Z Ho (/3)G Ho ^ H (P;-t-i/3) = Z H (/3)G H ^ Ho (P;t), (64) 

which, as expected, connects the generating function 
Gh () ^h of the work done by the "forward" quench to 
the one Gh^h of the "backward" quench, both originat- 
ing from a thermal state of the corresponding pre-quench 
hamiltonian (i.e., Hq and H respectively). Upon invert- 
ing the Fourier transform which defines the generating 
function, Eq. (64) yields 



Ph ->h(W) 
P h ^ Ho (-W) 



Zh(P) 
Zh W)' 



(65) 



which expresses the Crooks relation for the present case. 
Analogously, Eq. (64) at t = yields the Jarzyn- 
ski relation 3 "' (e'^) = e ~ p(FH ' FH o\ where F H = 
— lnZff(/3) represents the free energy of the system 
with hamiltonian H in equilibrium at an inverse temper- 
ature /3. 



VI. CONCLUSIONS 

In this paper we discussed the statistics of the work 
done by a quantum quench of the gap (mass) of a sys- 
tem of non-interacting bosons, which capture the low- 
energy properties of a variety of experimentally relevant 
cases. We derived the expressions for the work probabil- 
ity distribution, focusing on the universal features which 
emerge close to the edge singularity at low values of the 
work and upon approaching the critical point of the post- 
quench hamiltonian. We also studied thermal effects due 
to a non-zero initial temperature. 

For a system of non-interacting bosons, the calcula- 
tion presented here first of all confirms the applicabil- 
ity of the analytical continuation from imaginary to real 
time of the predictions derived in the context of bound- 
ary field theory : . In fact, an analogous analytic contin- 
uation has been shown to require some care in the case 
of the Ising chain in a transverse field quenched across 
its critical point ,J . By taking advantage of the relation- 
ship with the statistical properties of a classical gaussian 
field fluctuating in a film geometry in d + 1 space dimen- 
sions, we demonstrated that universal features emerge 
after a quantum quench of non-interacting bosons in d- 
dimensions upon approaching its critical (gapless) point, 
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which are in fact linked to the critical Casimir effect char- 
acterizing the corresponding classical system. In order 
for the features discussed here to emerge, the size L of 
the system has to be large compared to the correlation 
length £ of the post-quench hamiltonian and both of them 
have to be large compared to the microscopic length scale 
a (e.g., a possible lattice spacing) which characterizes the 
system. From the experimental point of view, measure- 
ments of the work statistics after a quantum quench can 
provide insight into the excitation spectrum of the quan- 
tum system under probe. Our analysis shows that in 
low-dimensional systems that can be described (cither 
exactly or approximately) by non-interacting, massive 
bosonic excitations, the work distribution is character- 
ized by a pattern of "fringes" which can be easily iden- 
tified. As a function of the work W, subsequent fringes 
are typically spaced by 2m, being m the post-quench gap, 
while their overall "amplitude" is set by the the square 
T of the fidelity (see Eq. (25)), which is exponentially 
suppressed as L increases. Correspondingly, the exten- 
sive mean value (W) of the work grows as oc L d , giving 
rise to a broad gaussian peak which carries most of the 
spectral weight and within which the various fringes 
are decreasingly visible (see Fig. 2). Accordingly, in or- 
der to detect the features we have highlighted here it 
is necessary to investigate systems which are large on a 
microscopic scale but still of finite size L ^> £ a. De- 
pending on the choice of the parameters of the quench m 
and mo, which determine the L-dependence of the overall 
amplitude T (see Eq. (25) and Fig. 5), it is possible to 
find an optimal combination of their values which clearly 
allows the detection of the edge singularity, as shown in 
Figs. 2 and 3. 

In addition, as the space dimensionality d of the model 
increases, the edge singularity associated to the low- 
est threshold at W = 2m becomes increasingly milder 
and, upon approaching the threshold from above, it ac- 
tually turns from a divergence for d < 2 to a vanish- 
ing distribution for d > 2 (see Fig. 3), with universal 
features encoded, inter alia, in the algebraic law which 
rules such an approach. In particular, while in one spa- 
tial dimension each peak is clearly distinguishable for 
the next one, as it decays to zero on a scale compa- 
rable to their separation, this is no longer the case in 
higher dimensions, where the weaker singularity results 
in a merging of the consecutive fringes. Note that, in con- 
trast to non-interacting fermionic systems where excited 
states cannot be more than singly or doubly occupied, 
in non-interacting bosonic systems multiple occupation 
of excited states leads to multiple excitation peaks. Ac- 
cordingly, the work done on the system can in princi- 
ple assume arbitrarily large values, though with a small 
probability, and therefore the support of the work prob- 
ability distribution is not bounded (see the discussion in 
Ref. 14). In the case of a quench occurring from an ini- 
tial thermal state, all the features highlighted above are 
lost and the probability distribution P(W) of the exten- 
sive work W takes the much less regular form reported in 



Fig. 7 for a certain choice of the parameters of the quench. 
Differently from the ordinary quench, there is no lower 
threshold for the work, which can assume arbitrarily large 
values, though with small probability. Accordingly, no 
universal features are expected to emerge at sufficiently 
small \W\. By comparing the mean work done during 
the quantum quench with the work that would be done 
in equilibrium during a reversible variation of the param- 
eter of the hamiltonian from its pre- to its post-quench 
value, one can study the associated irreversible entropy 
production AS^, which is expected to increase signif- 
icantly upon approaching the critical point 1 '. Rather 
generally it turns out that the pronounced maximum 
showed by ASi IT for narrow quenches and as a function 
of the quench parameter can be actually characterized in 
terms of the equilibrium specific-heat exponent and that 
a dimensional crossover can be observed upon varying 
the temperature of the initial state. An analysis along 
the one of Ref. 14 of the properties of large fluctuations 
around the broad gaussian distribution which develops 
upon increasing L is left for future investigations. 

The possibility to measure experimentally the work 
statistics would enhance our understanding of the excita- 
tion spectrum after a quantum quench and would provide 
an intriguing check of universality which is the corner- 
stone of the proposed connection between classical equi- 
librium physics in film geometries (and the associated 
critical Casimir effect) and non-equilibrium processes ex- 
emplified by the quantum quench problem. Even though 
no experimentally viable technique is presently available 
for measuring the statistic P(W) for many-body systems 
suitable extensions of recent proposals ' ' might provide 
experimental access to its generating function G(t). 
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Appendix A: Derivation of G(t) for a single 
harmonic oscillator. 

In this appendix, we consider a single harmonic oscil- 
lator described by the hamiltonian 

H=\p 2 + 1 -u>y, (Al) 
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in terms of the momentum p and the position q oper- 
ators and we calculate the characteristic function G(t) 
of the work done by a quantum quench of its frequency 
from u}q to w. The initial state is considered to be either 
the ground-state l^o) corresponding to frequency ujq or 
a thermal state at inverse temperature /3, always at fre- 
quency ujq. This problem is a special case of the one 
thoroughly addressed in Ref. 45 where the frequency is 
assumed to be any function of time w(r). Here, exploit- 
ing the simple step-like dependence of u> on the time in 
the quantum quench protocol, we will first present a cal- 
culation of G(t) for the case in which the initial state is 
the ground-state I'J'o) by expanding the latter in the basis 
of eigenstates of the post-quench hamiltonian H. Then 
we will repeat this calculation by using the coherent state 
formalism which has the advantage of being elegant, sig- 
nificantly simpler and straightforwardly generalizable to 
the thermal case, which we consider at the end of this 
appendix. 

Let us start with the calculation of G(t) by using the 
expansion of the initial state l^o) m Eq. (6) on the eigen- 
basis {|n)} of the post-quench hamiltonian: 



G(t) 



-(«+*) iw *|(n|¥o}| S 



(A2) 



where we used the fact that for a single harmonic oscilla- 
tor H\n) = (n + l/2)uj\n) and E — ojq/2. The first step 
is the calculation of the amplitudes (n|^ ): since both 
the pre- and post-quench hamiltonians are quadratic in 
the position and momentum operators q and p, the rela- 
tion between pre- and post-quench creation-annihilation 
operators ao, aj and a, a) (which are linear combinations 
of q and p) is linear. More specifically it is given by the 
single-boson version of the well-known Bogoliubov trans- 
formation 



a,Q = /ia + vw with \[i\ 



(A3) 



where the condition on fi and v is due to the requirement 
of having canonical commutation relations for the two 
pairs of creation-annihilation operators. The Bogoliubov 
coefficients /i, v for the quantum quench problem are real 
and given by 11 




(A4) 



The amplitudes (n^o) can be obtained readily by solv- 
ing the recurrence relation they satisfy. The latter can 
be derived by taking into account that aol^o) — and 



|n) = (at)" |0)/%M which yield: 

(n|*„) = ^<0|a n -Vao-w4)l*o) 
Vn! 



-[IV 



-[IV 



(n- 21* 



(OK/xa 1 + va)a T1 



(n-2|*o}-^ 2 (n|*o} ) 



(A5) 



where we made use of the relation Eq. (A3) between a 
and ao which implies the commutation relations [a, a ] = 
[x and [a n ,aQ] = n[ia n 
relation 



We thus find the recurrence 



(n|*n) = -A 



n- 1 



(n-2|* 



with 



_ V UJQ — LJ 

A = — = ■ — , 

fi u +0J 



(A6) 



(A7) 



which allows us to determine the sequence (n|*o) from 
its first two elements. A direct calculation based on 
the fact that the ground-state wave function *gs(<z) 01 
the harmonic oscillator with frequency us is v E'gs('?) 
exp(-wg 2 /2), yields (0|* ) = V^M and (1|* ) = 0. 
The latter relation implies that (n|*o) = whenever n 
is odd, which is expected because the wave function as- 
sociated with | Wo) is even under spatial parity g — > — q, 
while those associated with \n) have parity (—1)". The 
solution of the recurrence relation Eq. ( A6) is easily found 
in terms of double factorials as 



(n\*a) = 



(-A)"/ 2 



(n-1)! 



(A8) 



By using the property T(k + 1/2) = (2k - l)!! % /7r/2 fe 
and (2k)\\ = 2 k T(k + 1) (see, e.g., 5.4.2 in Ref. ), the 
solution for even n is 



(n|* > = (-A)" /2 



/Vrrr f + 1 



1/2 



(A9) 



while (n|*o) = for all odd n. Substituting into Eq. (A2) 
we have 



G(t) 



_ e*"° t/2 y r(n+i) , n 



'7T/X 



E 

71 = 



iujt \ 2n 



(A10) 



The sum of this series can be calculated by using the 
integral representation of T and the result is 



G(t) = 



e i(wo-w)i/2 
[iVl ~ A 2 e" 2i 



_ gi("o-")*/2 



1 - A 2 



^ _ ^2Q—2iuit ' 

(All) 
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We now repeat the derivation of Eq. (All) by us- 
ing coherent states and then we extend it to the ther- 
mal case. The first observation is that the initial state 
for the post-quench hamiltonian is a so-called squeezed 
vacuum, defined as a state of the form S(£)|0) where 
S(£) = exp (£*a 2 — £ a ^ 2 )] is the unitary squeezing op- 
erator and |0) is the vacuum state after the quench. This 
is a consequence of the Bogoliubov relations Eq. (A3) . In 
fact the squeezing operator S — 5(re 4 *) with coshr = fi 
and e^ sinh r — v performs exactly the Bogoliubov rota- 
tion according to a = SaS^, from which follows that 
the state 5|0) is annihilated by ao since ao(5|0}) = 
(SaSt)(S|0)) = Sa\0) = 0. Therefore |* ) = S^|0> or, 
equivalently, 



|*o> = (l-A 2 ) 1 / 4 



exp 



A 

— a 
2 



t2 



|0), 



(A12) 



as it can be shown for example by normal-ordering the 
operator S ( ). One can verify directly that ao given 
by Eq. (A3) annihilates the state in Eq. (A12) using the 
identity [a, /(at)] = /'( a t). 

The characteristic function G(t) can now be calculated 
by using known results from the theory of coherent and 
squeezed coherent states. Introducing in (\E , | e_4fft |^'o) 
of Eq. (6) a resolution of the identity in the basis of co- 
herent states of the post-quench hamiltonian, we get 



G(t) 



(*o|e~ 



iHt 



(A13) 



where \z) indicates a coherent state {a\z) = z\z) with 
complex z) and d 2 z = dRez dlmz. The action of the 
evolution operator on the coherent states is simply 



e- im \z) 



(A14) 



as one can easily verify by expanding the coherent state 
1 2) on the basis of the eigenstates of H. Accordingly, 
from Eq. (A13) one finds 



G(t) = e^ "")'/ 2 / — {^ \ze- iut ){z\^o}- (A15) 

7T 



The overlap between a coherent state \a) = D(a)\0) and 
a squeezed coherent state |/3;£) = S(£)D(j3)\Q), where 
D(a) = exp(aat — a*a) is the unitary displacement op- 
erator, is given by 0(1 



(z 

2 V 



2 Q -2iuit 



+ z* 2 ) 



a\z) — z\z) along with the standard relation (z|0) 
g-JI 2 ! determined by their normalization. 
Substituting Eq. (A17) in Eq. (A15) we find 

G{t) = / d 2 z exp 

TT/i J 

(A18) 

which is simply a double gaussian integral that can be 
evaluated straightforwardly, giving the final result (All). 

Note that the calculation presented above for 
( v I / o|e _lfft | v I / o) can also be done in imaginary time, i.e., 
for (^ r ol e_ff ' R |^ r o) which represents the partition function 
of the corresponding classical model in a film of width 
R with boundary states both equal to |\&o)j an d pro- 
vides the analytical continuation of the above result for 
t = —iR. The only difference is that, due to the standard 
normalization of coherent states, we now have 



-HR\ 



5 -io;H-||z| 2 (l-e 



\ze-" R ) (A19) 



instead of Eq. (A14). 

Apart from the overall phase factor e l ( LU o-u)t/2 
corresponds to the difference between the ground-state 
energies of the pre- and post-quench hamiltonian, the rest 
of Eq. (All) is a periodic function of t with period 2uj 
and as such it can be expanded as a Fourier series, from 
which wc deduce that the probability distribution Eq. (9) 
is given by an equidistant sequence of <5-peaks at positions 
W n = (u- Wo)/2 + 2nuj, n = 0,l, 2, 3, ... (Fig. 1). More 
precisely, by expanding Eq. (All) as a Fourier series in 
t (or Taylor expanding in powers of X 2 e~ 2lujt ) we recover 
Eq. (A10). 

We now consider a quantum quench from wo to u> but 
starting from a thermal initial state with inverse tem- 
perature j3, instead of the ground state. In this case the 
calculation in the Fock basis is significantly more difficult 
than before since knowledge of all transition amplitudes 
between eigenstates of the pre-quench hamiltonian with 
those of the post-quench one is required . The coherent 
state method instead remains relatively simple. 

The characteristic function G(t) is now given by 



G(t) 



Tr{c-P H °e iHot e- 



iHi 



} 



Tr {e-^o} 



(A20) 



By introducing a resolution of the identity in the coherent 
state basis of the post-quench hamiltonian, one finds 



(a\(3;0 = — e 
V A* 



^(l"| 2 + l/3| 2 )-3L( ra * 2 -^*/3 2 -2a*/3) 



(A16) 



where, as before, fi — coshp and v — e^sinhp. A 
squeezed vacuum is a squeezed coherent state with /3 = 
and so in our case, according to the above relations, we 
have 



(*l*0> 



1 



exp 



M 2 + Az* 2 



(A17) 



Alternatively, the last expression can be derived directly 
from Eq. (A12) and the definition of a coherent state 



Tr{e (_/3+i * )ff °e^ iff *} = J — 



, z U-0+U)Ho e -iHt, z s 



As in the ordinary case, the action of the evolution op- 
erator e~ lHt on \z) is simply given by Eq. (A14), but in 
order to proceed further we need to introduce once more 
a resolution of the identity, this time in the coherent state 
basis of the pre-quench hamiltonian 

Tr{e (_/3+ii)ffo e _jfft } = 



d 2 z d 2 C 



(z\e { - 0+u)Ho |0 (C\e-' lHt \z) , (A22) 
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where the states \Q are coherent states in the pre-quench 
basis (i.e., a \() = CIO) an d squeezed coherent states 
in the post-quench basis. The action of the operator 
c (-l3+it)H on j s now a ]^ g0 s i m piy given by Eqs. (A14) 

and (A19) and yields 



Tr{e (-,s+it)ffo e^ Ht } 

d 2 2d 2 C _l| C |2 a _ c -2^0 



e -Pu /2 e i(u -u)t/2 x 



-e 2 



\z\(e(-P +i V w °)(C\ze- iwt ), 
(A23) 



FIG. 8. Probability distribution P(W) for a thermal quan- 
tum quench in a single harmonic oscillator with cj — 0.1, 
Uq = 1, and ft = 0.5. This choice of tj and uiq is the same as 
the one of the quench reported is Fig. 1, which the present 
one has to be compared with. The red dots indicate peaks 
that are of purely thermal origin. 



which is a quadruple gaussian integral involving over- 
lap amplitudes between coherent and squeezed coherent 
states given by Eq. (A16). After some algebra and taking 
into account the normalization factor 1/Tr{e~^ }, we 
finally hnd 



G th (t) = 



e »(wo-w)t/ 2 (i _ e _ ^"°)\/l - A 2 



(A24) 



One can verify that this last expression agrees with the 
more general one of Ref. 45, Eq. (17). Note that in 
the limit /3 — > oo this expression renders Eq. (All), as 
expected. In contrast to the ordinary quench case in 
Eq. (All), the Fourier series of Eq. (A24) consists of dou- 
ble frequency oscillations since both e 1 " * and e iujt are in- 
volved. In order to highlight the emergence of novel peaks 
of thermal nature, we report in Fig. 8 the probability dis- 
tribution function P(W) obtained from Eq. (A24) for the 
single harmonic oscillator quenched from an initial state 
with ujq = 1 and ft — 0.5 to a final state with uj = 0.1, 
i.e., for the same values of parameters cuq and ui as those 
of the ordinary quantum quench reported in Fig. 1. Each 
vertical line indicate the occurrence of a (^-function at the 



I 

corresponding frequency, with the (integrated) amplitude 
reported in the vertical axis. In particular, the peaks 
highlighted in red are of purely thermal origin and are 
therefore absent in Fig. 1, together with all the peaks oc- 
curring for energies below the lower threshold (u — uiq) /2 
which characterize the ordinary quantum quench. In ad- 
dition, compared to the latter case (see Eq. (All)), e IW * 
in Eq. (A24) appears not only in even powers as before 
(reflecting the fact that the initial state \^q) contained 
excitations of the post-quench hamiltonian only in pairs) , 
but also in odd power terms (meaning that the thermal 
initial state contains all levels of excitations), with the 
latter decaying for small temperatures (/3 — > oo). 
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